%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Illustrate implications of counter-factuals for young workers
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

decomp=[];
decomp_l=[];

% XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
% SECOND-ORDER APPROXIMATION:
% XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

juebar=nanmean(jue_sim(:,:,1));
juibar=nanmean(jui_sim(:,:,1));
jeubar=nanmean(jeu_sim(:,:,1));
jeibar=nanmean(jei_sim(:,:,1));
jiebar=nanmean(jie_sim(:,:,1));
jiubar=nanmean(jiu_sim(:,:,1));

% Decomposition of U
%first-order terms
bue=-(jeubar.*(jiebar + jiubar).^2 + jeibar.*jiubar.*(jiebar + jiubar))./(jeubar.*jiebar + jeibar.*jiubar + jeubar.*jiubar + jiebar.*juebar + jiebar.*juibar + jiubar.*juebar).^2;
beu=(juebar.*(jiebar + jiubar).^2 + jiebar.*juibar.*(jiebar + jiubar))./(jeubar.*jiebar + jeibar.*jiubar + jeubar.*jiubar + jiebar.*juebar + jiebar.*juibar + jiubar.*juebar).^2;
bei=(juebar.*jiubar.^2 + jiebar.*(juebar + juibar).*jiubar)./(jeubar.*jiebar + jeibar.*jiubar + jeubar.*jiubar + jiebar.*juebar + jiebar.*juibar + jiubar.*juebar).^2;
bui=-(jeubar.*jiebar.*(jiebar + jiubar) + jeibar.*jiebar.*jiubar)./(jeubar.*jiebar + jeibar.*jiubar + jeubar.*jiubar + jiebar.*juebar + jiebar.*juibar + jiubar.*juebar).^2;
biu=(jiebar.*(jeibar.*juebar + juibar.*(jeibar + jeubar)))./(jeubar.*jiebar + jeibar.*jiubar + jeubar.*jiubar + jiebar.*juebar + jiebar.*juibar + jiubar.*juebar).^2;
bie=-(jiubar.*(jeibar.*juebar + juibar.*(jeibar + jeubar)))./(jeubar.*jiebar + jeibar.*jiubar + jeubar.*jiubar + jiebar.*juebar + jiebar.*juibar + jiubar.*juebar).^2;

% Second-order terms
bue2=(2.*(jiebar + jiubar).^3.*(jeubar - jeibar.*(jiebar./(jiebar + jiubar) - 1)))./((jiebar + jiubar).*(juebar + (jiebar.*juibar)./(jiebar + jiubar)) + (jiebar + jiubar).*(jeubar - jeibar.*(jiebar./(jiebar + jiubar) - 1))).^3;
beu2=-(2.*juebar.*(jiebar + jiubar).^3 + 2.*jiebar.*juibar.*(jiebar + jiubar).^2)./(jeubar.*jiebar + jeibar.*jiubar + jeubar.*jiubar + jiebar.*juebar + jiebar.*juibar + jiubar.*juebar).^3;
bei2=-(2.*juebar.*jiubar.^3 + 2.*jiebar.*(juebar + juibar).*jiubar.^2)./(jeubar.*jiebar + jeibar.*jiubar + jeubar.*jiubar + jiebar.*juebar + jiebar.*juibar + jiubar.*juebar).^3;
bui2=(2.*jeubar.*jiebar.^3 + 2.*jiubar.*(jeibar + jeubar).*jiebar.^2)./(jeubar.*jiebar + jeibar.*jiubar + jeubar.*jiubar + jiebar.*juebar + jiebar.*juibar + jiubar.*juebar).^3;
biu2=-(2.*jiebar.*(jeibar + jeubar + juebar).*(jeibar.*juebar + jeibar.*juibar + jeubar.*juibar))./(jeubar.*jiebar + jeibar.*jiubar + jeubar.*jiubar + jiebar.*juebar + jiebar.*juibar + jiubar.*juebar).^3;
bie2=(2.*jiubar.*(jeubar + juebar + juibar).*(jeibar.*juebar + jeibar.*juibar + jeubar.*juibar))./(jeubar.*jiebar + jeibar.*jiubar + jeubar.*jiubar + jiebar.*juebar + jiebar.*juibar + jiubar.*juebar).^3;

% Decomposition of LFPR
%first-order terms
bue_l=-((jeibar - juibar).*(jeubar.*jiebar + jeibar.*jiubar + jeubar.*jiubar))./(jeubar.*jiebar + jeibar.*jiubar + jeubar.*jiubar + jeibar.*juebar + jeibar.*juibar + jeubar.*juibar + jiebar.*juebar + jiebar.*juibar + jiubar.*juebar).^2;
beu_l=((jeibar - juibar).*(jiebar.*juebar + jiebar.*juibar + jiubar.*juebar))./(jeubar.*jiebar + jeibar.*jiubar + jeubar.*jiubar + jeibar.*juebar + jeibar.*juibar + jeubar.*juibar + jiebar.*juebar + jiebar.*juibar + jiubar.*juebar).^2;
bui_l=-((jeibar + jeubar + juebar).*(jeubar.*jiebar + jeibar.*jiubar + jeubar.*jiubar))./(jeubar.*jiebar + jeibar.*jiubar + jeubar.*jiubar + jeibar.*juebar + jeibar.*juibar + jeubar.*juibar + jiebar.*juebar + jiebar.*juibar + jiubar.*juebar).^2;
bei_l=-((jeubar + juebar + juibar).*(jiebar.*juebar + jiebar.*juibar + jiubar.*juebar))./(jeubar.*jiebar + jeibar.*jiubar + jeubar.*jiubar + jeibar.*juebar + jeibar.*juibar + jeubar.*juibar + jiebar.*juebar + jiebar.*juibar + jiubar.*juebar).^2;
biu_l=((jeibar + jeubar + juebar).*(jeibar.*juebar + jeibar.*juibar + jeubar.*juibar))./(jeubar.*jiebar + jeibar.*jiubar + jeubar.*jiubar + jeibar.*juebar + jeibar.*juibar + jeubar.*juibar + jiebar.*juebar + jiebar.*juibar + jiubar.*juebar).^2;
bie_l=((jeubar + juebar + juibar).*(jeibar.*juebar + jeibar.*juibar + jeubar.*juibar))./(jeubar.*jiebar + jeibar.*jiubar + jeubar.*jiubar + jeibar.*juebar + jeibar.*juibar + jeubar.*juibar + jiebar.*juebar + jiebar.*juibar + jiubar.*juebar).^2;

% Second-order terms
bue2_l=(2.*(jeibar - juibar).*(jeibar + jiebar + jiubar).*(jeubar.*jiebar + jeibar.*jiubar + jeubar.*jiubar))./(jeubar.*jiebar + jeibar.*jiubar + jeubar.*jiubar + jeibar.*juebar + jeibar.*juibar + jeubar.*juibar + jiebar.*juebar + jiebar.*juibar + jiubar.*juebar).^3;
beu2_l=-(2.*(jeibar - juibar).*(jiebar + jiubar + juibar).*(jiebar.*juebar + jiebar.*juibar + jiubar.*juebar))./(jeubar.*jiebar + jeibar.*jiubar + jeubar.*jiubar + jeibar.*juebar + jeibar.*juibar + jeubar.*juibar + jiebar.*juebar + jiebar.*juibar + jiubar.*juebar).^3;
bei2_l=(2.*(jeubar + juebar + juibar).*(jiubar + juebar + juibar).*(jiebar.*juebar + jiebar.*juibar + jiubar.*juebar))./(jeubar.*jiebar + jeibar.*jiubar + jeubar.*jiubar + jeibar.*juebar + jeibar.*juibar + jeubar.*juibar + jiebar.*juebar + jiebar.*juibar + jiubar.*juebar).^3;
bui2_l=(2.*(jeibar + jeubar + jiebar).*(jeibar + jeubar + juebar).*(jeubar.*jiebar + jeibar.*jiubar + jeubar.*jiubar))./(jeubar.*jiebar + jeibar.*jiubar + jeubar.*jiubar + jeibar.*juebar + jeibar.*juibar + jeubar.*juibar + jiebar.*juebar + jiebar.*juibar + jiubar.*juebar).^3;
biu2_l=-(2.*(jeibar + jeubar + juebar).^2.*(jeibar.*juebar + jeibar.*juibar + jeubar.*juibar))./(jeubar.*jiebar + jeibar.*jiubar + jeubar.*jiubar + jeibar.*juebar + jeibar.*juibar + jeubar.*juibar + jiebar.*juebar + jiebar.*juibar + jiubar.*juebar).^3;
bie2_l=-(2.*(jeubar + juebar + juibar).^2.*(jeibar.*juebar + jeibar.*juibar + jeubar.*juibar))./(jeubar.*jiebar + jeibar.*jiubar + jeubar.*jiubar + jeibar.*juebar + jeibar.*juibar + jeubar.*juibar + jiebar.*juebar + jiebar.*juibar + jiubar.*juebar).^3;


% Construct actual and counterfactual U rates:
jei_alt=jei;jui_alt=jui;jue_alt=jue;jeu_alt=jeu;jie_alt=jie;jiu_alt=jiu;
%the original are stored under jAB_0

%focus on young only:
wbar=zeros(1,11);
wbar(1)=.5;
wbar(6)=.5;

for kk=1:2
    if kk==1
    jei=jei0;jui=jui0;jue=jue0;jeu=jeu0;jie=jie0;jiu=jiu0;
    elseif kk==2
    jei=jei_alt;jui=jui_alt;jue=jue_alt;jeu=jeu_alt;jie=jie_alt;jiu=jiu_alt;
        
    end
    
    juebar=mean(jue);
    jeubar=mean(jeu);
    juibar=mean(jui);
    jeibar=mean(jei);
    jiebar=mean(jie);
    jiubar=mean(jiu);
    
    %GET THE COMPONENTS OF UNEMPLOYEMENT TO SECOND-ORDER
    
    %%%%%%%%%%% U %%%%%%%%%%%
    eps_ue=kron(bue,ones(length(jue),1)).*(jue-kron(juebar,ones(length(jue),1)));
    eps_ue2=1/2*kron(bue2,ones(length(jue),1)).*(jue-kron(juebar,ones(length(jue),1))).^2;
    epsw_ue=kron(wbar,ones(length(jue),1)).*(eps_ue+eps_ue2);%+(epsep_ue+epset_ue+epseq_ue+epsue_ui+epsue_ei+epsue_iue+epsue_ilf));
    
    eps_eu=kron(beu,ones(length(jue),1)).*(jeu-kron(jeubar,ones(length(jue),1)));
    eps_eu2=1/2*kron(beu2,ones(length(jue),1)).*(jeu-kron(jeubar,ones(length(jue),1))).^2;
    epsw_eu=kron(wbar,ones(length(jue),1)).*(eps_eu+eps_eu2);%+(epsep_ue+epsep_ui+epsep_ei+epsep_iue+epsep_ilf));
    % NOTE: we ignore the cross-order terms between the separation E-X transitions rates
    % they are too small to matter
    
    eps_ei=kron(bei,ones(length(jue),1)).*(jei-kron(jeibar,ones(length(jue),1)));
    eps_ei2=1/2*kron(bei2,ones(length(jue),1)).*(jei-kron(jeibar,ones(length(jue),1))).^2;
    epsw_ei=kron(wbar,ones(length(jue),1)).*(eps_ei+eps_ei2);%+(epsep_ei+epset_ei+epseq_ei+epsue_ei+epsui_ei+epsei_iue+epsei_ilf));
    
    eps_ui=kron(bui,ones(length(jue),1)).*(jui-kron(juibar,ones(length(jue),1)));
    eps_ui2=1/2*kron(bui2,ones(length(jue),1)).*(jui-kron(juibar,ones(length(jue),1))).^2;
    epsw_ui=kron(wbar,ones(length(jue),1)).*(eps_ui+eps_ui2);%+(epsep_ui+epset_ui+epseq_ui+epsue_ui+epsui_ei+epsui_iue+epsui_ilf));
    
    eps_ie=kron(bie,ones(length(jue),1)).*(jie-kron(jiebar,ones(length(jue),1)));
    eps_ie2=1/2*kron(bie2,ones(length(jue),1)).*(jie-kron(jiebar,ones(length(jue),1))).^2;
    epsw_ie=kron(wbar,ones(length(jue),1)).*(eps_ie+eps_ie2);%+(epsep_ie+epset_ie+epseq_ie+epsue_ie+epsui_ie+epsei_ie+epsei_ilf));
    
    eps_iu=kron(biu,ones(length(jue),1)).*(jiu-kron(jiubar,ones(length(jue),1)));
    eps_iu2=1/2*kron(biu2,ones(length(jue),1)).*(jiu-kron(jiubar,ones(length(jue),1))).^2;
    epsw_iu=kron(wbar,ones(length(jue),1)).*(eps_iu+eps_iu2);%+(epsep_iu+epset_iu+epseq_iu+epsue_iu+epsei_iu+epsui_iu+epsiue_iu));
    
    
    %%%%%%%%%%% LFPR %%%%%%%%%%%
    epsl_ue=kron(bue_l,ones(length(jue),1)).*(jue-kron(juebar,ones(length(jue),1)));
    epsl_ue2=1/2*kron(bue2_l,ones(length(jue),1)).*(jue-kron(juebar,ones(length(jue),1))).^2;
    epswl_ue=kron(wbar,ones(length(jue),1)).*(epsl_ue+epsl_ue2);%+(epsep_ue+epset_ue+epseq_ue+epsue_ui+epsue_ei+epsue_iue+epsue_ilf));
    
    epsl_eu=kron(beu_l,ones(length(jue),1)).*(jeu-kron(jeubar,ones(length(jue),1)));
    epsl_eu2=1/2*kron(beu2_l,ones(length(jue),1)).*(jeu-kron(jeubar,ones(length(jue),1))).^2;
    epswl_eu=kron(wbar,ones(length(jue),1)).*(epsl_eu+epsl_eu2);%+(epsep_ue+epsep_ui+epsep_ei+epsep_iue+epsep_ilf));
    % NOTE: We ignore the cross-order terms between the separation E-X transitions rates
    % they are too small to matter
    
    epsl_ei=kron(bei_l,ones(length(jue),1)).*(jei-kron(jeibar,ones(length(jue),1)));
    epsl_ei2=1/2*kron(bei2_l,ones(length(jue),1)).*(jei-kron(jeibar,ones(length(jue),1))).^2;
    epswl_ei=kron(wbar,ones(length(jue),1)).*(epsl_ei+epsl_ei2);%+(epsep_ei+epset_ei+epseq_ei+epsue_ei+epsui_ei+epsei_iue+epsei_ilf));
    
    epsl_ui=kron(bui_l,ones(length(jue),1)).*(jui-kron(juibar,ones(length(jue),1)));
    epsl_ui2=1/2*kron(bui2_l,ones(length(jue),1)).*(jui-kron(juibar,ones(length(jue),1))).^2;
    epswl_ui=kron(wbar,ones(length(jue),1)).*(epsl_ui+epsl_ui2);%+(epsep_ui+epset_ui+epseq_ui+epsue_ui+epsui_ei+epsui_iue+epsui_ilf));
    
    epsl_ie=kron(bie_l,ones(length(jue),1)).*(jie-kron(jiebar,ones(length(jue),1)));
    epsl_ie2=1/2*kron(bie2_l,ones(length(jue),1)).*(jie-kron(jiebar,ones(length(jue),1))).^2;
    epswl_ie=kron(wbar,ones(length(jue),1)).*(epsl_ie+epsl_ie2);%+(epsep_ie+epset_ie+epseq_ie+epsue_ie+epsui_ie+epsei_ie+epsei_ilf));
    
    epsl_iu=kron(biu_l,ones(length(jue),1)).*(jiu-kron(jiubar,ones(length(jue),1)));
    epsl_iu2=1/2*kron(biu2_l,ones(length(jue),1)).*(jiu-kron(jiubar,ones(length(jue),1))).^2;
    epswl_iu=kron(wbar,ones(length(jue),1)).*(epsl_iu+epsl_iu2);%+(epsep_iu+epset_iu+epseq_iu+epsue_iu+epsei_iu+epsui_iu+epsiue_iu));
       
    % XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
    
    % ordering of decomposition:  UE, EU, EI, UI, IU, IE
    decomp{kk}=[sum(epsw_ue,2) sum(epsw_eu,2) sum(epsw_ei,2) sum(epsw_ui,2)  sum(epsw_iu,2) sum(epsw_ie,2) ];
    decomp_l{kk}=[sum(epswl_ue,2) sum(epswl_eu,2) sum(epswl_ei,2) sum(epswl_ui,2)  sum(epswl_iu,2) sum(epswl_ie,2) ];
end

% contribution of changes in LF shares
epsw_lf=kron(Ubar,ones(length(jue),1)).*(w-kron(wbar,ones(length(jue),1)));

lil=kron(mean(Lr{1})./mean(Lr_agg{1}),ones(length(jue),1));

%%%%%%%%%%%%%%%%%%%% FIGURES %%%%%%%%%%%%%%%%%%%%%%%%%

figure, 
subplot(311), 
bar(1:length(nber),-1+max(0,5*nber),1,'EdgeColor',[.7,.7,.7],'FaceColor',[.7,.7,.7],'LineStyle','none','Clipping','on'),xlim([1 length(nber)])
h1=gca;
set(h1,'YAxisLocation','right','Color','white','XTickLabel',[]);
set(h1,'Fontsize',8)
h2=axes('Position',get(h1,'Position'));
set(h2,'Fontsize',8)

tt=(100*[sum(decomp_l{1},2)+mean(wbar*Lr{1}')]);
tt2=100*(sum(decomp_l{2},2)+mean(wbar*Lr{1}'));

plot(hpfilter(tt,10),'LineWidth',1);
hold on, plot(hpfilter(tt2,10),'LineWidth',2,'Color',[1 0 0],'LineStyle','-');
set(h2,'YAxisLocation','left','Color','none')
ylim(h1,[2 3]), set(h1,'Yticklabel',[]), set(h1,'YTick',[],'XTick',[])
ylabel('ppt')
ylim(h2,[50 70])
set(gca','Xticklabel',num2str((1980:10:2015)'))
set(gca,'Xtick',12*4+1:12*10:length(tt)), xlim([1 length(tt)]),
set(h2,'XLim',get(h1,'XLim'),'Layer','top')
legend('boxoff')
h1 = legend('$LFP$','$\widehat{LFP}$','Location','Best');
set(h1,'Interpreter','latex','FontName','Arial','FontSize',7)
legend('boxoff'), title('Young')

subplot(312)
bar(1:length(nber),-1+max(0,5*nber),1,'EdgeColor',[.7,.7,.7],'FaceColor',[.7,.7,.7],'LineStyle','none','Clipping','on'),xlim([1 length(nber)])
h1=gca;
set(h1,'YAxisLocation','right','Color','white','XTickLabel',[]);
set(h1,'Fontsize',8)
h2=axes('Position',get(h1,'Position'));
set(h2,'Fontsize',8)
plot(100*hpfilter(sum(decomp{1},2)+mean(wbar*Ur{1}'),10),'LineWidth',1);
hold on, plot(100*hpfilter(sum(decomp{2},2)++mean(wbar*Ur{1}'),10),'LineWidth',2,'Color',[1 0 0],'LineStyle','-');
set(h2,'YAxisLocation','left','Color','none')
ylim(h1,[2 3]), set(h1,'Yticklabel',[]), set(h1,'YTick',[],'XTick',[])
ylabel('ppt')
ylim(h2,[8 20])
set(gca','Xticklabel',num2str((1980:10:2015)'))
set(gca,'Xtick',12*4+1:12*10:length(tt)), xlim([1 length(tt)]),
set(h2,'XLim',get(h1,'XLim'),'Layer','top')
legend('boxoff')
h1 = legend('$U$','$\widehat{U}$','Location','Best');
set(h1,'Interpreter','latex','FontName','Arial','FontSize',7)
legend('boxoff')

diff=100*[sum(decomp{1},2)-sum(decomp{2},2)]';
subplot(313), plot(hpfilter(diff,10^6),'Linewidth',2)
hold on, plot(0*diff,'--black')
ylabel('ppt','FontSize',8), ylim([-1.5 1.75])
set(gca,'FontSize',8)
h3 = legend('$U-\widehat{U}$','Location','NorthEast');
set(h3,'Interpreter','latex','FontName','Arial','FontSize',7)
set(gca','Xticklabel',num2str((1980:10:2015)'))
set(gca,'Xtick',12*4+1:12*10:length(tt)), xlim([1 length(tt)]),
legend('boxoff')

%Data for online appendix
trend_yg=hpfilter(diff,10^6);

